Finite Differences#
강좌: 수치해석 프로젝트
Taylor Expansion#
Taylor series를 이용하면 함수를 쉽게 근사화 할 수 있다.
여기서 \(T.E\) 는 Truncation error 이다. 만약 1차 미분까지만으로 근사화 한 경우 이 오차는 \((\Delta x)^2, (\Delta x)^3, ...\) 과 같이 간격 \(\Delta x\) 의 고차 항으로 구성되어 있다.
\(\Delta x\) 가 작은 경우 오차는 Leading error 항과 그보다 매우 작은 항의 합이다.
수치 미분#
Taylor expansion 을 이용하면 수치 미분을 구할 수 있다.
Forward Difference
Backward Differnce
Central Difference
다음 두 식을 빼보자.
그 결과는 다음과 같다.
2차 미분
위의 두 식을 더한 후 \(2 f(x_j)\) 를 빼보자.
일반적인 방법
여러 Taylor expansion의 합차를 이용해서 원하는 미분항의 근사식을 구한다.
예제#
\(f(x)=\sin(x)\) 에 대해서 수치 미분 결과를 비교하자.
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
plt.style.use('ggplot')
plt.rcParams['figure.dpi'] = 150
def forward_diff(f, x, dx):
# Forward difference
return (f(x+dx) - f(x)) / dx
def backward_diff(f, x, dx):
# Backward difference
return (f(x) - f(x-dx)) / dx
def central_diff(f, x, dx):
# Central difference
return (f(x+dx) - f(x-dx)) / (2*dx)
def compute(dx):
x = np.linspace(0, 2*np.pi, 101)
f = np.sin
# Compute first derivatives
exact = np.cos(x)
fd = np.array([forward_diff(f, xi, dx) for xi in x])
bd = np.array([backward_diff(f, xi, dx) for xi in x])
cd = np.array([central_diff(f, xi, dx) for xi in x])
return x, exact, fd, bd, cd
def plot(dx):
x, exact, fd, bd, cd = compute(dx)
# Plot
plt.plot(x , exact)
plt.plot(x, fd)
plt.plot(x, bd)
plt.plot(x, cd)
plt.legend(['Exact', 'Forward Difference', 'Backward Difference', 'Central Difference'])
plt.xlabel(r'x')
plt.ylabel(r"$f'(x)$")
plt.title("Comparison of finite difference @ $\Delta x$={} $\pi$".format(dx/np.pi))
plot(0.2*np.pi)

plot(0.1*np.pi)

정확도 비교#
Forward / Backward difference 는 1차 정확도 (\(O(\Delta x)\))
Central difference 는 2차 정확도(\(O((\Delta x)^2)\))
def error(dx):
_, exact, fd, bd, cd = compute(dx)
err_fd = np.linalg.norm(fd - exact)
err_bd = np.linalg.norm(bd - exact)
err_cd = np.linalg.norm(cd - exact)
return err_fd, err_bd, err_cd
# Change delta x
dxs = [10**(-n) for n in range(1, 6)]
err_fd, err_bd, err_cd = [], [], []
for dx in dxs:
fd, bd, cd = error(dx)
err_fd.append(fd)
err_bd.append(bd)
err_cd.append(cd)
# Plot error (log-log scale)
plt.loglog(dxs, err_fd, marker='o')
plt.loglog(dxs, err_bd, marker='x')
plt.loglog(dxs, err_cd, marker='d')
plt.legend(['Forward Difference', 'Backward Difference', 'Central Difference'])
plt.xlabel(r'$\Delta x$')
plt.ylabel('Error')
Text(0, 0.5, 'Error')

Modified wavenumber#
정확도를 평가하는 또 다른 방법
임의의 함수는 Fourier expansion을 이용하면 Sinusoidal 함수로 표현 가능
\(n\) 이 커질수록 주기가 짧아짐 : High frequency
High frequency 항을 모사하기 위해서는 \(\Delta x\) 를 줄이거나 정확도가 높아야 함
Finite difference 의 정확도에 따라 Low frequeny 항은 잘 모사할 수 있으나 High frequency 항은 왜곡될 수 있음
x = np.linspace(0, 2, 101)
legends = []
for i in range(1, 4):
plt.plot(x, np.sin(np.pi*i*x))
legends.append("Mode {}".format(i))
plt.legend(legends)
<matplotlib.legend.Legend at 0x7fe0a6b842e0>

Modified wavenumber
간단한 Harmonic 함수에 대해 분석
\([0, L]\) 영역을 N개의 격자로 나누었다고 생각하자. 이때 waveumber \(k\)를 다음과 같이 고려하자
이론적인 미분은
Finite difference
격자점은
- Central differnce에 대해 Harmonic 함수 적용
즉
kh = np.linspace(0, np.pi, 101)
mkh = np.sin(kh)
plt.plot(kh, kh, linestyle='--')
plt.plot(kh, mkh)
plt.xlabel('k h')
plt.ylabel("k'h")
Text(0, 0.5, "k'h")

plot(0.5*np.pi)

plot(np.pi)
